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ABSTRACT 

Context. The observed relationship between mass, age and rotation in open clusters shows the progressive development 
of a slow-rotator sequence among stars possessing a radiative interior and a convective envelope during their pre-main 
sequence and main-sequence evolution. After 0.6 Gyr, most cluster members of this type have settled on this sequence. 
Aims. The observed clustering on this sequence suggests that it corresponds to some equilibrium or asymptotic condition 
that still lacks a complete theoretical interpretation, and which is crucial to our understanding of the stellar angular 
momentum evolution. 

Methods. We couple a rotational evolution model, which takes internal differential rotation into account, with classical 
and new proposals for the wind braking law, and fit models to the data using a Monte Garlo Markov Chain (MCMC) 
method tailored to the problem at hand. We explore to what extent these models are able to reproduce the mass and 
time dependence of the stellar rotational evolution on the slow-rotator sequence. 

Results. The description of the evolution of the slow-rotator sequence requires taking the transfer of angular momentum 
from the radiative core to the convective envelope into account. We find that, in the mass range 0.85-1.10 Mq, the core¬ 
envelope coupling timescale for stars in the slow-rotator sequence scales as Quasi-solid body rotation is achieved 

only after 1-2 Gyr, depending on stellar mass, which implies that observing small deviations from the Skumanich law 
(P oc \/i) would require period data of older open clusters than is available to date. The observed evolution in the 0.1-2.5 
Gyr age range and in the 0.85-1.10 Mq mass range is best reproduced by assuming an empirical mass dependence of the 
wind angular momentum loss proportional to the convective turnover timescale and to the stellar moment of inertia. 
Period isochrones based on our MCMC fit provide a tool for inferring stellar ages of solar-like main-sequence stars 
from their mass and rotation period that is largely independent of the wind braking model adopted. These effectively 
represent gyro-chronology relationships that take the physics of the two-zone model for the stellar angular momentum 
evolution Into account. 

Key words. Stars: rotation - Stars: evolution - Stars: late-type - open clusters and associations: general - open clusters 
and associations: individual: Pleiades, M34, M35, M37, Hyades, Praesepe, Coma Berenices NGC6811 NGC6819 


1. Introduction 


In the past decade, stellar rotation has received a great 
deal of attention, both because o f its promising potential 
as a precise stel lar age estimator ( Barnes|2007 2010 Mei- 


bom et al. 12015 1 and the large amount of good quality data 


increasingly available, either as a byproduct of planetary 
transit searches or from d edicated stellar rotation observa- 
tional programs (see, e.g. Herbst fc Mundt]|2005[ Hodgkin 


et al. 


et al. 


2006 


2008 


][ 200 ^ 

von Braun et al.||2005[ |Messina[|2007[ ‘ 1M essina| 


principle be 


2010||. Rotation periods of young stars can in 
derived with a precision of 10“^ 


using photo¬ 
metric monitoring, with surface differential rotation being 
probably the most important limiting factor for individual 
old stars (see, e.g. Epstein fc Pinsonneault|[2014 l. 

The realization by Barnes (20031 of the existence of a 
well-defined sequence of slowly rotating stars in the young 
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open clusters color-period diagram has provided a criterion 
for establishing an empirical relationship between mass, age 
and rotation, which can be used to infer the age of stars 
that a re sufficiently old to belong to that sequence. |Barnes| 
(20031 effectively disentangled the convergence of stellar 


periods towards a unique slow-rotator sequence, seen as 
a progressive reduction of the scatter that is due to the 
presence of fast rotators (from very young clusters to the 
Hyades), from the rotational evolution of stars in this se¬ 
quence, which dominates the col or-period diagr am from the 
Hyades onwards. Subsequently, Barnes (20071 provided a 
calibration of his rotation-mass-age relationship using the 
observed rotation-mass distribution in open clusters and 
the solar age, rotation period and {B — V). More recently, 
Barnes & Kim| (20101 provided evidence of a connection 
between the empirical rotational period functional depen¬ 
dence on (B — V) and the con vective turnover timescale 
T. Based on this, jBarnes (20101 derived a nonlinear model 
formulated in terms of the Rossby number (Ro = P/t) and 
two dimensionless empirical constants, which describes the 
rotational evolution towards the slow-rotator sequence. 
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In this work we focus on the rotational evolution of stars 
already on the slow-rotator sequence and study how well 
the current models are able to reproduce such evolution. 
We focus on the age interval between 0.1 and 2.5 Gyr, for 
which data with the richest details is available, providing 
strong constraints on the models. The observed clustering 
on this sequence in the mass-age-period space suggests that 
it corresponds to some equilibrium or asymptotic condition 
that still lacks a complete theoretical interpretation and 
may be crucial to our understanding of the stellar angular 
momentum evolution. 

Data that has become available in recent years has 
revealed an evolutionary behavior of the slow-rotator se¬ 
quence that s eems difficult to reconcile with the original 
Barnes (20031) idea of a factorization of the time and mass 

_ TNotably, in their stud y of rotation in M35 and 

M34, Meibom et al. ( 2009[ 2011b[ | found that while G-type 
stars on the slow-rotator sequence spin down with a rate 
close to the Skumanich P oc ^/i law, K-type stars spin down 
more slowly. 

On the theoretical side, recent work concentrated on 
improving the description of the magnetic braking due to 
coupling between the stellar wind and magnetic fiel d s gen- 
erated by an internal stellar dynamo. Matt et al. (20121 
performed two-dimensional axisymmetric magnet ohy dro- 
dynamic simulations to compute steady-state solutions for 
solar-like stellar winds from rotating stars with dipolar 
magnetic fields, from which they derived a semi-analytic 
formulat i on for the external torque on the star. |Gallet fc] 
Bouvier ( 2013|) derive d the angular momentum loss rate 
using the Matt et al. (20121 equation for the Alfven ra¬ 
dius together with a dynamo prescription that relates the 
stellar magnetic field to stellar rotation, as well as a wind 
prescription that relates the mass-loss rate to the stellar an¬ 
gular velocity, obtained from current theory and calibrated 
to the present-day Sun. Given the current uncertainties on 
the s tellar magnet i c field s and on the wind mass-loss, how¬ 
ever, Matt et al.l (20151 derived a scaling for the depen- 


2. Data 

For this work we use literature data of stars in open clus¬ 
ters that show a well-defined slow-rotator sequence, easily 
identifiable by eye as an over-density in the period-color 
diagram, also known in t he literature as the “I-sequence” 
(Barnes 20031 or “ridge” (Kovacs 20151. The age at which 
this sequence begins to form is still rather undetermined, 
but it is evident that this is already well defined at the age 
of the Pleiades and present in all older open clusters. A 
summary is reported in Tablej^ The data set includes clus¬ 
ters from 0.1 to 2.5 Gyr for which [B — V) measurements 
in the relevant range are readily available. The age interval 
is sampled in an optimal way, in the sense that data for 
open clusters not considered here would not improve the 
age sampling in a significant way. Since our method relies 
on an empirical fitting of the slow-rotator sequence to start 
with (see Sect.j^, data for field stars are not considered. 
Young loose stmlar associations are also not considered in 
this work. 

Rotational pe riods for the Pleiades are taken from |Hart^ 
man et al. (2010). This data set includes 241 single stars in 
common with the comp i lation of probable Pleiades mem¬ 
bers by Stauffer et al.| (20071, which provides near- and 
mid-infrared photometry together with a compilation of ho- 
mogenise d Johnson B and V m agnitudes. For the Pleiades, 
following Stauffer et al.| (19981 we adopt an age of 120 Myr. 


Meibom et al. ( 2009|) provide rotational periods and UB- 
VRl photometry for 310 member s of M35. The age o f M35 
has been estimated to 150Myr (von Hippel et al. 2002), 


175 Myr dBarrad o y Navascu6s e t al. 200] j), an d 180 Myr 
(Kalirai et al. 2003|). h'ollowing Meibom et al. ( 2009| , we 
adopt an age of 150 Myr. 

Rotational periods and reddening correct ed {B — V)rj 
for M34 are taken from Meibom et al. (201 lb I. The sample 
contains periods for 118 stars. From the rotational period of 
slow-rotator sequence stars, adopting the Skumanich law, 
Meibom et al. (2011b| estimated an age of 240 Myr, which 
we also adopt here. 


dence of the stellar wind braking on Rossby number and 
an empirically-derived scaling with stellar mass and radius. 

In this work, we include these various wind braking laws 
in the classical two-zone description of the stellar rotational 
evolution and investigate to what extent they reproduce 
the slow-rotator sequence evolution. The models are fit¬ 
ted to the data by means of a Monte Garlo Markov Ghain 
(MGMG) method tailored to the case at hand. Using the pa¬ 
rameters derived by the MGMG model fitting to the data, 
we derive period isochrones and tracks from 0.1 to 5 Gyr 
that are largely independent of the specific wind braking 
law adopted. These effectively represent gyro-chronology 
relationships that take the effects of differential rotation 
in the period evolution into account. 

The paper is organized as follows. In Sect .1^ we discuss 
the data used in the present analysis. In Sect.j^we present 
an empirical description of the slow-rotator sequence, which 
also illustrates the current state-of-the-art, together with a 
novel non-parametric empirical fit of its evolution. In Sect.|^ 
we present the MGMG fitting method of two-zone rota¬ 
tional evolution models with different wind braking pre¬ 
scriptions. In Sect, ^the results of the MGMG fitting to the 
data are presented and discussed. The conclusions are in 

Sect.m 


Hartman et al. (20091 provide a clean sample of rota- 
tional periods for 372 stars in M37. Photometry for this 


cluster is fr om Kalirai et al. (20011. Following Hartman 
et al. (20081 we adopt an age of 550 Myr. 


(20111 


oers 


of 


The rotational period data set of Delorme et al. 
includes 56 members of the Hyades and 35 mem 
Praesepe. They also provide BV photometry collected from 
the lit erature. Following the discussion in [Delorme et al.| 
(20111, we adopt an age of 600Myr for both clusters. 


Collier Cameron et al. (20091 present a photometric sur¬ 
vey oTTotationTates'lirthe^Coma Berenices cluster con¬ 
taining data for 27 members with available BV photome¬ 
try. Following their discussion, we adopt an age for Coma 
Berenices equal to the Hyades and Praesepe age (600 Myr) 
within the uncertainties. 


Meibom et al. (201 la I provide rotational periods for 71 
cluster members of NGC6811. Photometry for this cluster 


is provided by panes et al.j (|2013|) . Following [Meibom et al. 
( [2011a and references therein) we adopt an age of 1 Gyr. 
Rotatio nal periods and phot ometry for NGC6819 is 


taken from Meibom et al. 


(20151, who derived a gyro age 
of 2.5 Gyr in agreement with the estimate by [Jeffries et al.[ 

dM^. 


It is worth noting that the cluster ages adopted here 
do not differ significantly from the isochrone ages adopted 
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Table 1. Clusters data used in this work. 


Name Age (Myr) references 


Pleiades 

120 



Hartman et al. 

( 

2010) 

M35 

150 



iVieibom et al. 

( 

2009) 

M34 

220 



Vieibom et al. 

1 

2TJTTb) 

M37 

550 



Hartman et a 


1 

2009) 

Hyades 

600 



Uelorme et al. 

( 

2011) 

Praesepe 

600 



Uelorme et al. 

( 

2011) 

Goma Berenices 

600 

Gollier 

Gameron et al. 

( 

2009) 

NGG6811 

1000 

Meibom et aT7 

(20TT 

31 

Janes et al. 

( 

2013) 

NGG6819 

2500 



Meibom et al. 

( 

mS) 


by Kovacs (20151. In the present work, isochrone and rota- 
tional ages are assumed to coincide. 


3. Empirical description of the slow-rotator 
sequence evolution 

3.1. Empirical models 

The main assumption in Barnes ( 2003[ 20071 is that ro¬ 
tation periods on the slow-rotator sequence can be repre¬ 
sented empirically by a relationship of the form 


P = g{t)f{B-V), 


( 1 ) 


where t is the stellar age and the colour B — V is taken as a 
proxy of the stellar mass. Other authors fitt ed the function 
P = P{t, B — V) to different data s ets (e.g. Meibom et al. 
2009 Mamajek fc HillenbrandS2008|), mainta ining t he same 
form and factorisation as in |Barnes| (| 2007| . Both iBarnesI 


(20081 derived their 
studied clusters ( e.g . 


kc \PoJ 


( 2 ) 


-•'-"-I? 


(3) 


which predicts that, asymptotically, P —)■ g(t )f(B—V) with 
g{t) = y/t and f{B — V) = The 


Barnes 


(20101 


The color-period diagrams for the clusters studied here 
are shown in Fig.j^ The ste llar mass is estimated using 


the color-mass relationship of Barnes & Kim (20101, based 
on the effective temper ature-color transformation of |Leje-| 


19981. Periods are scaled to the square 


(20071 and Mamajek & Hillenbrair 
color-period relationship on a few wel 
a Pers ei, Pleiades, Hyades) while the fit by Meibom et al. 
( 20091 is based o n a single cluster (M35). 

Barnes (2003|) proposed that the mass dependence of the 
rotation periods in open clusters can be simply attributed 
to the moments of inertia of either the whole star or of 
the outer convection zone. This was later found ina dequate 
by Barnes & Kiml (20101, and led Barnes (2010) to the 
formulation of a dinerent non-linear empirical model, which 
provides the gyro-age of a star as 


where kc and fc/ are two dimensionless constants, con¬ 
strained by the observations, and Pq is the initial rotation 
period, taken as the period on the zero-age main sequence. 
In the slow-rotator sequence limit we are concerned with, 
this semi-empirical formula gives 


une et al.| ( |1997[ 
root of age (in Myr), so that data points would overlap in 
diagrams at different ages if the rotational evolution were 
exactly Skumanich-type. The slow-rotator sequence is eas¬ 
ily identifiable in all clusters, with the only exception of 
the Pleia des a nd M35, for which some ambiguity may arise 
(see Sect. 3.2). From the comparison of the observed slow- 
rotator sequence at different ages it is evident that: 

— Stars with M < 1.0 Mq slow down at a slower rate than 
predicted by the Skumanich law until t ~ 0.55 Gyr; 

— From 0.55 to 2.5 Gyr, stars slow down at a faster rate 
than predicted by the Skumanich law, except for the 
NGG6811 stars in the 0.8 < M/Mq < 0.9 mass range; 

— There is no evidence of stars of mass M < 0.8 Mq on the 
slow-rotator sequence at ages younger than 0.24 Gyr. 

The current observed evolution of the slow-rotator se¬ 
quence from 0.1 to 2.5 Gyr shows that the factorization ex¬ 
pressed by Eq. 0 does not hold, at least before 0.55 Myr, 
because such a relationship would imply that the shape of 
the slow-rotator sequence is maintained throughout its evo¬ 
lution. 

The data available contains periods for stars on the 
slow-rotator sequence with mass as low as M « 0.5Mq 
for ages «0.6 Gyr, but no data for stars with M < O.SMq 
is available at later ages. 

In Fig.we also report the f(B — V) empirical func¬ 
tions obtained by Barnes|(2007l, Meibom et al.|(2009l, and 
Mamajek fc Hillenbrand ] 

of Eq. (I 2 I) from |Barnes 


as well as the inversion 
(20101) and its asymptotic limit in 


Eq. 0.^ ?he comparison with observations shows that such 
relationships can describe the slow-rotator sequence only 
on limited age and mass ranges. 

With the exception of stars of 0.8 < M < 0.9 Mq in 
NGG6811, the slow-rotator sequence lies above the « 0.55 
Gyr sequence in the P/Vi vs. {B — V) d iagram. In the 
framework of the two-zone model (Sect. 4.1), this behavior 


model gives a description of the mass-dependent evolution 
of fast-rotating stars to the slow-rotator sequence. An alter- 
n ative descriptio n is the metastable dynamo model (MDM) 
of Brown (2014), which assumes a stochastic transition be¬ 
tween two diff erent dynamo regimes and is rooted in an 
original idea of Barnes (2003). 


could be attributed to the transport of angular momentum 
from the stellar core to the envelope in the earlier evolution 
after the ZAMS. We shall verify this in Sect.|^ 

3.2. Non-parametric empirical fitting 

In order to overcome the limitations of the empirical model¬ 
ing discussed in the Sect. |3.l] we perform a non-parametric 
fit to the slow-rotator sequence. To this end, we need a 
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Pleiades ( 120 Myr) 

Moss 

1.1 1.0 0.9 0.8 0.7 0.6 



M35 ( 150 Myr) 

Moss 

1.1 1.0 0.9 0.8 0.7 0.6 



M34 ( 240 Myr) 

Moss 

1.1 1.0 0.9 0.8 0.7 0.6 



M37 ( 550 Myr) 

Moss 

1.1 1.0 0.9 0.8 0.7 0.6 



Hyades - Preaesepe - Coma Berenices ( 600 Myr) 

Moss 

1.1 1.0 0.9 0.8 0.7 0.6 



NGC681 1 (1000 Myr) 

Moss 

1.1 1.0 0.9 0.8 0.7 0.6 



NGC6819 (2500 Myr) 

Moss 

1.1 1.0 0.9 0.8 0.7 0.6 



(B-V)„ 


Fig. 1. Rotation period scaled to the square root of age (in Myr) vs. {B — V). Periods of stars defining the slow-rotator sequence 
are represented with black dots, all others with grey dots. Our non-parametric ht of the slow-rotator sequence is shown as blue 
thick solid lines, with the Pleiades and M37 fits repeated in all pan els as referenc e (black thin solid line and grey thick solid line 
respectively). The dotted line represent the f{B — V) function of |Barnes| (|2007|l (B07), the dashed line that of|Meibom et al.| 
(2009) (M09), the dot-dashed that of Mamajek & Hillenbrand ([2008| (M08) , the dash-triple-dotted line the inverted Barnes (|2010|) 
relationship (Eq.[^ BIO), the long-dashed line the [Barnes r pM r asymptotic slow-rotator sequence limit (Eq.|^ BlOa). 
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Table 2. Rotational periods from the non-parametric fit. a is the standard deviation, p the Pearson test probability that the 
residual distribution is normal. 


cluster 

age 

1.10 

1.05 

1.00 

0.95 

0.90 

M/Mq 

0.85 

0.80 

0.75 

0.70 

0.65 

0.60 

(7 

V 

Pleiades 

120 

2.09 

2.75 

3.53 

4.42 

5.39 

6.41 

7.42 





0.35 

0.63 

M35 

150 

3.12 

3.62 

4.16 

4.82 

5.64 

6.46 






0.76 

0.80 

M34 

240 

3.18 

3.92 

4.72 

5.60 

6.57 

7.65 

8.65 

9.48 

10.11 

10.62 


0.53 

0.82 

M37 

550 

5.21 

5.91 

6.50 

7.08 

7.70 

8.42 

9.32 

10.32 

11.19 

11.91 

12.43 

0.94 

0.68 

Hyades 

600 



7.62 

8.35 

9.07 

9.82 

10.64 

11.57 

12.26 

12.72 

13.04 

0.51 

0.52 

Praesepe 

600 


6.53 

7.36 

8.16 

8.90 

9.55 

10.02 

10.48 

11.09 

11.77 


0.63 

0.89 

ComaBerenices 

600 

6.68 

7.49 

8.12 

8.70 

9.29 

10.07 

10.75 

11.07 

11.23 

11.37 


0.52 

0.26 

Hydes-Praesepe-Coma 

600 

6.12 

6.95 

7.65 

8.31 

8.99 

9.68 

10.38 

11.20 

11.77 

12.22 

12.58 

0.73 

0.68 

NGC6811 

1000 

7.07 

8.58 

9.90 

10.93 

11.30 

11.03 

10.64 





0.46 

0.45 

NGC6819 

2500 

15.46 

17.41 

19.12 

20.75 

21.87 

22.07 






0.88 

0.92 


Pleiades 


U35 


M34 





M37 



AP 


Hyades • Praesepe • Coma Berenices 



AP 


NGC6811 



NGC6819 



Fig. 2. Comparison of the empirical distribution function of the non-parametric fit residuals with the normal distribution function. 
Residuals are normalized to the standard deviation. 


criterion for selecting the stars belonging to this sequence. 
For the older clusters, the selection is quite simple as only 
some outliers and a few remaining fast-rotators need to 
be excluded. For the younger clusters, on the other hand, 
the slow-rotator sequence over-density is still easily identi¬ 
fiable in the color-period diagram, but the separation with 


stars approaching the slow-rotator sequence is somewhat 
arbitrary and prone to subjective choices. Slow-rotator se¬ 
quence m embership ma y become even more ill-defined if, 
following Barnes (2010), we consider the sequence as an 
asymptotic limit. 
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We note, however, that non-parametric fits to the older 
slow-rotator sequences produce residuals whose distribution 
can be approximated with a normal distribution. To trace 
back the sequence to the younger clusters with a criterion 
as homogeneous as possible, we set the width of the slow- 
rotator sequence as the maximum width that produces a 
distribution of residuals with a normality probability (see 
below) of at least 30% (see Fig.j^and Ta blef^. We use the 
local polynomial regression fit^OESS, Cleveland et al. 


1992[) with smoothing para meter between 0.8 and 1 . 0, and 


a Pearson normality test (Moore 1986 Thode Jr. 20021, 
as implemented in the statistical package R. In order to 
improve the fitting at the higher and lower end of the se¬ 
quence, where the curvature is higher and a non-parametric 
fit would require a much denser data set, we use Eq. ([^ as a 
normalization function, adjusting kj in order to mimic the 
slow-rotator sequence at « 0.6 Gyr. 

Our non-parametric fits are shown in Fig.j^ where the 
stars selected as members of the slow-rotator sequence are 
also highlighted. The non-parametric fit is used to derive 
periods for a grid of stellar masses, reported in Table|^ 
together with the standard deviation and the normality 
test probability. The comparisons of the empirical distri¬ 
bution functions of the residuals with normal distributions 
are shown in Fig.[^ 

Our selection of stars belonging to the slow-rotator se¬ 
quence in younger clust ers is more re strictive than that 
adopted, for instance, by Barnes (20101. It should be con¬ 
sidered, however, that the aim of the present work is to 
study in detail the rotational evolution of the slow-rotator 
sequence, under the assumption that this represents some 
equilibrium state or asymptotic behavior, deliberately ig¬ 
noring fast rotators or stars still approaching the slow- 
rotator sequence. We also ignore periods much larger than 
the periods of the slow-rotator sequence, assuming that 
they are either incorrect or that these stars are not clus¬ 
ter members. In other words, we assume that the peak of 
the over-density of slow rotators in the color-period diagram 
corresponds to the maximum probability of having reached 
the equilibrium/asymptotic state and that stars in or near 
such a state have periods distributed normally around the 
peak. Such an approximation has the advantage of allow¬ 
ing a characterization of the width of the sequence using 
the standard deviation, which can then be used as input 
in our parametric fit described in Sect.|^ It turns out that 
such an approximation describes remarkably well the dis¬ 
tribution in the older clusters, and that, by an appropriate 
selection, it can be extended to the younger clusters as well. 
For our purposes, it can safely be assumed that the stan¬ 
dard deviation of the sequence includes both observational 
uncertainties and the intrinsic physical width due to the 
ongoing convergence onto the slow-rotator sequence. 


4. Two-zone models Monte Carlo Markov-chain 
fitting to the data 

4.1. Two-zone rotational evolution models 


In our modeling of the rotational evolution of solar-like stars 
(whose structure is characterized by an inner radiative re¬ 
gion, surrounded by a convective envelope), we adopt the 
theoretical framework of two-zone models (TZMs; see, e.g. 
MacGregor fc Brenner|[T991[ [Keppens et al.|[l995| |Allain| 


1998 Spada et al.|201l| . The main assumption of the TZM 


Article number, page 6 of 


16 


is that the radiative interior and the convective envelope of 
the star rotate as rigid bodies, with angular velocities ilc 
and Oe, respectively; the rotational state of the whole star 
at a given time is thus completely specified by these two 
quantities. 

During the pre-main sequence, the rotational evolution 
is dominated by evolutionary changes to the stellar struc¬ 
ture (i.e., overall contraction, appearance and growth of a 
radiative core), and the interaction with the circumstellar 
disc. From the zero age main sequence onwards, the interior 
structure changes very little, and the rotational evolution 
essentially depends on the interplay between the angular 
momentum loss at the surface, due to the magnetized stel¬ 
lar wind, and the angular momentum exchange between 
core and envelope. These main physical ingredients are in¬ 
cluded in our models as follows: 


The radiative core grows at the expenses of the sur¬ 
rounding en velope, thus s ubtracting from it the angular 
momentum (Allain|l998l: '^McRc^e, where Me and Rc 
are the mass and the radius of the core, respectively (in 
the following, we adopt the dot notation for derivatives 
with respect to time); 

The star-disc interaction is assumed to be very efficient, 
capable of keeping the surface angular velocity constant 
during the entire disc lifeti me, raise (the so-called disc¬ 
locking approximation, see Koenigl|[l991 1; 

The rate of angular momentum loss from the surface 
due to the magnetized stellar win d is specified by the 
wind braking law Jwb (see Section 
Over the coupling timescale Tep, 

lar momentum AJ = — ° (0^ — Og) is transferred 

T Ie 

fro m the radiative interior to th e convective envelope 
(see MacGregor fc Brenner|[T9^ |. 


4.2 


for details); 


the amount of angu- 


The TZM equations for and Ge are therefore: 
- For t < Tdisc: 


/A =+2/3Mci?^ Oe - icG,; 

G. =0; 


(4) 


- For t > Tdisc: 


/cGc= +2/3McRl^e 
Ietle= -2l3MeRl^e 


^J/'Tcp ic^c 

T At//"Tep /gOe 


Jwb- 

(5) 


In the equations above, we have introduced the moments of 
inertia, R and R, of the radiative zone and of the convec¬ 
tive envelope, respectively. Note the appearance of terms 
involving the time derivatives of the moments of inertia 
(/gGe, etc.), as required by angular momentum conserva¬ 
tion when significant structural changes take place. 

We begin the integration of equations (|^ and © when 
the star is still fully convective, when it is reasonable to 
assume that the star is rotating as a solid body. We thus 
set the initial conditions in terms of a single parameter, the 
initial period Pq: 


Gc(to) = Ge(<o) 


27r 
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We extract the evolution of stellar structure quantities, 
Rc, Me, etc., from stellar models of appropriate mass, cal- 
cnlated nsing the Yale Rotational stellar Evolntion Code 


(YREC), i n its non-rotational configuration (Demarque 
et al. ]2008| . The stellar models were calculated with initial 


composition and mixing le ngth parameter oivilt cali brated 
on the Sun. Assnming the Grevesse & Sauval (1998 1 value 
of the solar metallicity, (i^/A) 0 ^GS 98 ~ 0.0230, our solar 
calibration gives Xq — 0.7039, Yq = 0.2773, Zq = 0.0188, 
C^MLT = 1.826. 

4.2. Wind braking laws 

In our TZM framework, the wind braking law prescribes 
the dependence of J^b on the stellar mass (either directly 
or indirectly, i.e., through a dependence on other stellar 
structure variables), and on the snrface angular velocity 
Dg. In this work, we foens on the following wind braking 
laws: 


1. The Kawaler (19881 law as modified by Chaboyer et al. 
([T995fl 


-K„ 


Jvjh = < 


\M,/Mq) 


-K„ 


R. 


/Rq \ 


1/2 


M,/Mq J 


n for Dg ^ flgat 


( 6 ) 


Note that the dependence on reproduces the empir¬ 
ical Skumanich law if the star rotates as a solid body. 
Our scaling of Kyj reqnires that the numerical value 
given in Sect. [^mnst be mnltiplied by 1.11 • 10 "*^ g cm^ s 
to convert it into cgs units. 

2. A hybrid braki ng law, which c ombines the classical 
dependence of Kawaler (|1988|) wit h the mass depen¬ 
dence snggested by Rarnes fc Kim (20101 and Barnes 


(20101 for the slow-rotator seqnence: 


Jwh 




(7) 


where /» and Iq are the moment of inertia of the whole 
star and of the Sun, respectively, and is au ad- 
jnstable parameter with the same scaling as above. 


3. The braking law proposed by Gallet & Bonvier (20131. 


Their starting point is the very general expression Jwb oc 
OeMwind i’a) where Mwind is the mass loss dne to the 
stellar wind and ta is the Alfven radius, i.e., the radius 
of the (approximately) spherical region around the star 


where the wind is d ominated by the magneti c field (We- 
ber fc Davis||1967| . [Gallet fc Bonvier (20131 use an 


ex- 


pression deri ved from the numerical simulation of [Matt | 


et al. (20121 to eliminate the Alfven radins from the 


braking law; the ir final resnlt is (see Sect. 3.3 of Gallet 
fc Bonvi^|2013 for details): 


Jwb = Kf fleMwind 


f?2 p2 


-\ 2 m 


Mwind + niRl 


Ri 


In the equation above, = 2GM*/i?* is the escape 
velocity at the stellar snrface. Bp is the snrface mag¬ 
netic field and Mwind is the stellar mass loss rate due 
to the wind. In the equation above, the constants Ki, 
K 2 , an d m are fixed accor ding to the nnmerical simula¬ 
tions of Matt et al. (20121; their valnes are: Ki = 1.30, 


K 2 = 0.0506, m = 0.2177. The values of B^, and Mwind 


are calcnla ted nsing the BOREAS snbroutine (Granmer fc 
Saar|20lT I with the filling factor slightly modified in or¬ 


der to rep roduce the average filling factor of the present 
Sun as in Gallet fc Bouvier ( 2013[). 


4. The braking law proposed by Matt et al. (20151. They 
derived a physically motivated scaling for the depen¬ 
dence of the wind braking law on Rossby nnmber and 
adopted an empirical scalin g with stellar mass and ra¬ 
dins. For our purposes, the Matt et al. (20151 braking 


law can be recast in the form: 


Jwh — 


with 






for Ro > Roq /X 
for Ro < Roq /X 


3.1 


Mq 


0.5 


and 


Ro, 


X = 


'© 


i?Og 


^sat'^ 

^qTq 


(9) 


( 10 ) 


( 11 ) 


Ddel 


Following Matt et al.| (120151, we adopt p = 2, which 


implies that this model also reproduces the empirical 


Skumanich (1972| law in the solid-body rotation limit. 
The satnTalion regime is equivalent to that introduced 


by Krishnamurthi et al. (1997); here we adopt x = 10 
as in Matt et al. (2015p (see also Sect.[^. 


4.3. MCMC fitting of the TZM parameters 

According to the TZM, the rotational evolution is com¬ 
pletely determined by the stellar mass, which also selects 
the appropriate backgronnd stellar model, and by five pa¬ 
rameters: the initial period Pq, the disc lifetime Tdiso the 
conpling timescale Tep, and the two parameters contained in 
the wind braking laws ([^ or 0 - We nse a MGMG approach 
to constrain the values of these parameters. 

The MGMG method is a powerful technique to obtain 
quantitative inferences on a set of model paramete rs p from 
the observa tional data d (see, e.g. Appendix A of Tegmark 
et al.|2004 for a concise introdnetion). In the case at hand, 
the vector of model parameters for the TZM with, e.g, the 
braking law ([^ is p = (Pq, "Tdisc, "Tep, flgat), while for 
each mass listed in Table the vector of data points con¬ 
tains the periods from the non-parametric fits at the varions 
cluster ages. For the models with M* = 1.00M©, we also 
consider the additional constraint based on the enrrent so¬ 
lar rotation rat^ fle(t 0 ) = 1.0 ± 0.1 XIq. 


(8) ^ We adopt the values t© = 4.57 Gyr, SI© = 2.903 • 10 ® s 
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Formally, Bayes’ Theorem relates the probability of 
the model, given the data, i.e., the posterior probabilit'^ 
7^(p|d), to the compatibility of the data with the model, 
i.e., the likelihood 7^(d|p), and the intrinsic likeliness of the 
model (prior probability, or simply prior), ’P(p): 

^(P|d) = P) ^(P)- (12) 

In the following, we ignore the normalization 7^(d), which 
is a constant with respect to the models, and we use the 
notation £(d; p) for the likelihood. In essence, the MCMC 
method consists of sampling the posterior probability by 
means of a quasi-random walk in the space of the parame¬ 
ters p, constructed as a chain of jumps guided by the values 
of the likelihood £(d;p). Each step in the chain depends 
only on the previous one (Markov Chain), and is partly 
stochastic (Monte Carlo). The steps are performed accord¬ 
ing to the following two-stage rule: 


1 . 


2 . 


Propose stage: Poid Ptriai = Poid + Ap, with Ap 
extracted from the jump function /(Ap); 
Acceptance/rejection stage: the probability to accept the 
jump is given by the Metropol is-Hastings rule (Metropo- 
lis et al.|l953 Hastings||1970| |: 


P accept — min 


l^(d; Ptriai) ^ 
/:(d;poid) ’ 


(13) 


The Metropolis-Hastings rule ensures that a trial step in¬ 
creasing the likelihood is always accepted (£(d; ptriai) > 
£(d;poid) => T’accept = 1), while at the same time even 
steps which reduce the likelihood retain a (proportionally 
small) non-zero probability to be accepted. This procedure 
is completely specified by the choice of the functions /(Ap) 
and £(d; p). 

We use a multi-dimensional Gaussian form for the jump 
function: 

/(Ap) oc , (14) 

3 


where £j are the jump sizes and the index j (up to 5 in 
our case) runs over the elements of the parameter vector 
p. The jump function is critical to ensure that the steps 
in the chain are optimally correlated, providing a useful 
sampling of the posterior probability: the optimal regime 
is in between 100% step acceptance (i.e., a purely random 
walk) and 100% step rejection. 

We define the likelihood in terms of the chi-square: 

£(d; m) = (15) 


The chi-square contains contributions from each cluster for 
which an estimate of the rotation period on the slow-rotator 
sequence at the mass considered is available: 


= E 


2=1 




(16) 


where £lrn,i = ^e(ti) is the angular velocity of the convec¬ 
tive envelope calculated from a run of the TZM with the set 

^ The standard notation V(A\B) denotes the conditional prob¬ 
ability of event A, given that event B is observed. 
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of parameters p, while are derived from Table 

as follows: 

O — — o 

^^ 0,2 - 7 - ^^ 0,2 ■ 

2 

In practical applications, running the MCMC procedure 
at the optimal acceptance regime (^ 20% of attempted 
steps accepted), is often hindered by the existence of cor¬ 
relations in the way the parameters influence the model. In 
our case, for example, we can anticipate the existence of de¬ 
generate pairs of values of Pq and ruigc, i.e., a longer initial 
period and a shorter disc interaction, and vice versa, com¬ 
pensating each other to give a very similar evolution at suffi¬ 
ciently late times (say, a few hundred Myr onwards). Apart 
from computational efficiency, strong correlations can lead 
to oversampling of the degenerate regions of the parame¬ 
ter space, resulting in very biased, and practically useless, 
chains . We d eal with this issue as recommended by |Tegmark| 
et al. (20041, by introducing a transformation from p to 
an uncorrelated vector u, and quantitatively checking the 
quality of the chain: 


— To transform to uncorrelated variables, we calculate the 
correlation matrix of the parameters 

C = (PP*) - (P) (p‘) 

over a portion of the chain (typically, every 100 steps), 
then diagonalize it to obtain the matrices A and R: 

C = RAR*; 


since C is by definition real and symmetric, the eigenval¬ 
ues are real and the matrix R is orthogonal. The linear 
transformation 


u = A-i/2Rt [p- (p)] 

results in the very useful properties: (u) = 0, (uu*) = 1, 
i.e., the variables Ui are (linearly) uncorrelated and nor¬ 
malized to one. Operationally, the vector p is trans- 
form ed i nto u before substitution into the jump func¬ 
tion (14l; the trial step is performed in the variables u, 
and trim transformed back if the step is accepted. Since 
u is normalized to one, the choice of the jump sizes is 
particularly simple: £j = £ k, 1. The correlation matrix 
C is recalculated, and the matrices A and R are updated, 
every ^100 steps. 

To evaluate the quality of the chain, we calculate the 
autocorrelation of each parameter pj , 

, . _ {pj{i -b fc)pj-(i)) - (pj-(f))^ 

^ kv]ii))-{V3id)Y ’ 


where Pj{i), etc., is the value of pj at the step i of the 
chain. Typically, Cj{k) decreases with k (its value at A: = 
1 being unity by definition): we define the correlation 
length of the chain as the value kj for which Cj{k) drops 
below 0.5 for the first time. The effective length of the 
chain Nj is thus the number of steps N divided by kj: 

N, = N/k,. 

This is a measure of the number of independent points 
in the chain: if Nj is sufficiently large, the average of 
Pj over the chain, along with its standard deviation, 
can be considered representative of the best-fitting value 
of Pj and its uncertainty, respectively. To ensure that 
we only use a portion of the chain sufficiently close to 
convergence, the first few hundred steps (burn-in phase) 
are discarded. 
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model KA model KS 




model KK 



model KB 



Fig. 3. Correlation of the parameter A™ for models KA, KS, KK, KB, MB, and Ki for model GB with mass. The scaling factor of 
Km is 1.11 • lO^^gcm^ s. The best linear fit is reported together with the uncertainties on the slope parameter. The linear Pearson 
correlation coefficient is given is each panel. 


5. Results and discussion 

The empirical description of Sect. [^implies that, once a star 
settles on the slow-rotator sequence, it f orgets most of its 
earlier rotational history. As suggested by |Barnes (20031 or 
Brown (20141, fast rotators may have a different magnetic 
configuration than slow rotators and migrate to the slow- 
rotator sequence after a change in their magnetic configu¬ 
ration. This scen ario, which still lacks clear observational 
support (but see Garraffo et al. 20151, implies that stars 


may settle on the slow-rotator sequence either after an ini¬ 
tial rotational evolution like those described in Sect. 14.l| and 
|4.2[ with parameters constant in time, or after an entirely 
different one, in which the param eters change with surface 
magnetic field configuration (e.g. Brown|[2014 |. In the first 
case, convergence on the slow-rotator sequence would be 
just a consequence of the dependence of the wind brak¬ 
ing law, including the effects of the saturation regime, which 
is implicitly assumed in building our starting conditions. In 
the second case, both the overall scale factor and the expo- 
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Table 3. Slow-rotator sequence two-zone Monte Carlo Markov-chain fitting results (see text for details on the different models). 
For each mass and model this table reports the resulting the fitted for models KA, KS, KK, KB, and MB or for model 
GB, and the fitted core-envelope coupling timescale Tcp. For both K and Tcp the table also gives the uncertainty and the effective 
chain length. The scaling factor of is 1.11 • lO'^^ gcm^ s. 


Mass {Mq) 

Model 


K 

ctk 

Nk 

Tcp (Myr) 

(Myr) 

Nr 

Tcp 


KA 

2.00 

6.03 

0.48 

12250 

85.5 

11.1 

8166 


KS 

2.03 

6.16 

0.48 

12250 

86.8 

10.7 

8166 

0.85 

KK 

1.72 

5.47 

0.40 

12250 

79.7 

10.6 

8166 

KB 

1.99 

4.78 

0.37 

12250 

84.6 

10.8 

8166 


MB 

2.12 

5.32 

0.41 

8166 

83.1 

11.5 

8166 


GB 

1.54 

2.42 

0.08 

12250 

85.0 

9.1 

8166 


KA 

1.21 

5.24 

0.50 

8166 

51.3 

9.4 

8166 


KS 

1.25 

5.32 

0.50 

8166 

52.6 

9.4 

8166 

0.90 

KK 

1.28 

5.41 

0.49 

8166 

53.5 

9.1 

8166 

KB 

1.61 

5.03 

0.48 

8166 

59.6 

9.3 

8166 


MB 

1.44 

5.28 

0.48 

8166 

52.3 

9.9 

8166 


GB 

1.09 

2.43 

0.10 

12250 

67.4 

7.1 

12250 


KA 

1.05 

4.90 

0.61 

8166 

34.9 

9.4 

8166 


KS 

1.11 

4.96 

0.62 

8166 

35.5 

9.1 

8166 

0.95 

KK 

1.45 

5.13 

0.52 

8166 

38.3 

7.9 

8166 

KB 

1.07 

4.32 

0.54 

8166 

35.1 

9.1 

8166 


MB 

1.11 

5.38 

0.65 

8166 

34.9 

9.4 

8166 


GB 

1.14 

2.16 

0.11 

8166 

49.0 

6.4 

8166 


KA 

1.81 

3.82 

0.64 

8166 

19.1 

8.5 

8166 


KS 

1.75 

3.76 

0.58 

8166 

19.0 

8.1 

8166 

1.00 

KK 

1.75 

3.89 

0.54 

8166 

20.0 

7.4 

8166 

KB 

1.62 

3.60 

0.58 

8166 

18.7 

8.1 

8166 


MB 

1.30 

4.51 

0.68 

8166 

17.3 

8.3 

8166 


GB 

2.11 

1.91 

0.13 

8166 

34.1 

5.7 

8166 


KA 

1.87 

3.79 

0.91 

8166 

15.6 

7.9 

8166 


KS 

1.90 

3.69 

0.77 

8166 

15.6 

7.5 

8166 

1.05 

KK 

1.81 

3.76 

0.87 

8166 

15.4 

7.8 

8166 

KB 

1.73 

4.69 

1.07 

8166 

15.3 

7.4 

8166 


MB 

1.71 

5.52 

1.23 

8166 

15.0 

8.3 

8166 


GB 

2.37 

1.85 

0.18 

8166 

23.0 

5.2 

8166 


KA 

3.45 

3.17 

0.91 

6125 

12.0 

6.6 

6125 


KS 

3.52 

3.40 

0.87 

6125 

12.8 

6.4 

6125 

1.10 

KK 

3.61 

3.65 

1.17 

6125 

12.9 

7.3 

6125 

KB 

3.25 

4.55 

1.43 

6125 

12.5 

6.9 

6125 


MB 

3.29 

6.97 

2.27 

6125 

12.4 

7.9 

6125 


GB 

2.29 

1.87 

0.28 

8166 

14.1 

4.3 

8166 


nent of 17 in the wind braking law may change at some stage 
of the evolution, which would require a different approach 
than that adopted here. In both cases, however, once on 
the slow-rotator sequence, the subsequent rotation period 
evolution is manifestly independent of the previous history, 
although the abundance of li ght elements may still depend 
on it (see, e.g. Bouvier|[2008 1. 


Our fitting, therefore, cannot constrain the stellar ro¬ 
tational history before the settling on the slow-rotator se¬ 
quence. In particular, it cannot constrain Pq ^-nd Tdisc even 
if the earlier rotational history follows one of the mod¬ 
els listed in Sect. |4.l| and |4.2| with parameters constant in 
time. Nevertheless, such parameters are necessary in our 


Article number, page 10 of 


modeling for advancing the rotational evolution from the 
birth-line to the slow-rotator sequence. We therefore treat 
Pq and Tdisc as nuisance parameters with reasonable priors 
and marginalize them out. The prior on Pq is unifor mly dis- 
tributed over the observed period range in the ONC ([Rebull 


et al. 2004). The prior on Tdisc is assumed to be an expo- 
nential d ecay with an e-folding time of « 5 Myr (see, e.g. 
Fig. 11 in Hernandez et al.||2008 l. 


In our modeling, Pq and Tdisc are degenerate with re¬ 
spect to the slow-rotator sequence evolution, in the sense 
that long Pq and short Tdisc can lead to the same period at 
the age of the Pleiades as short Pq and long Tdisc- Therefore, 
after checking compatibility by letting all parameters vary. 
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Fig. 4. Comparison of the mass dependence of models KA (solid 
black line), KB (red dashed line), and MB (blue dot-dashed line) 
at t = 0.55 Gyr. All quantities are in solar units. 


we set Tdisc = 3Myr and marginalize Pq out for each mass 
and model. It is worth stressing that this is just a convenient 
way of setting the initial conditions for the slow-rotator se¬ 
quence evolution and that no inference can be made on the 
actual distribution of Pq and ruisc- 

The duration of the phase of the early stellar evolution 
that is affected by the saturation regime (Eq.|^ or de¬ 
pends on the assumptions made on the saturation thresh¬ 
old. As a consequence, the star may leave the saturation 
regime before or significantly after its settling on the slow- 
rotator sequ ence. In fact, if we conside r the relationship 
proposed by Krishnamurthi et al. (19971: 


flgat — flgat©; 
TQ 


(17) 


with risat© = 10 fl©, or, equivalently, t he saturat ion r egime 
in the formulation of Matt et al. (20151 (see Sect. 4.21, stars 
of approximately solar mass will leave the saturation regime 
before settling on the slow-rotator sequence (starting from 
the Pleiades age in our modeling), but stars of lower mass 
will stay in the saturation regime for a significant part of 
their earlier evolution on the slow-rotator sequence. For this 
reasons, we adopt the strategy of testing different assump¬ 
tions on risat j rather than trying to fit it as a free parameter. 
We therefore consider Eq. (|y with: 


— no saturation (model KA); 

^ Usat = ion© for all masses (model KS); 

— risat as in Eq. ( [l7| (model KK). 

Note that an implicit mass dependence (through r) enters 
both model KK, as a ©-dependent saturation threshold, and 
Eq.g (hereafter model KB), as a factor in the wind angu¬ 
lar momentum loss. On the other hand, Eq. ^ (hereafter 
model MB) contains a complex mass dependence that is 
both explicit and implicit, through the dependence on r 
and i?*, and the ©-dependence of the saturation threshold, 
respectively. 

The wind braking law (|^ (model GB hereafter) con¬ 
tains, in principle, no adjustable parameters. However, in 
order to compare it with models derived from Eqs. @,0, 



Fig. 5. ©cp vs. M for models KA (orange circle), KS (red pen¬ 
tagon), KK (plum diamond), MB (purple hexagon), KB (green 
square) and GB (blue triangle). Results from different models 
at the same mass are slightly shifted in abscissa for readabil¬ 
ity. The solid line represents the relation ©cp oc obtained 

from model KB with an average (K-w) = 4.5. 


and (0, where is treated as a free parameter, we as¬ 
sume iCi variable and check whether the fitted Ki shows 
some residual correlation with mass. 

By fitting or Ki, we both test how well the models 
describe the mass dependence of the rotational evolution, 
and correct the mass dependence using the observational 
data. 

After setting the initial conditions, we are left with two 
free parameters, p = (©cp, for the models KA, KS, KK, 
KB, and MB and p = (©cp, Ki) for model GB. We fit these 
parameters using the MGMG method described in Sect. 4.3 
and the data of Table |2] From this data set we exclude 
data for stellar mass < 0.85 M©, since, for our purpos^ 
the sampling in age in this mass range is still insuflicienlM 
We also exclude the NGG6811 periods, although we verify 
a posteriori (see below) the consistency of our results with 
the NGG6811 data. The exclusion of the NGG6811 periods 
is due to the irregular behavior in the empirical description 
of the evolution of the slow-rotator sequence of this cluster 
in comparison to the general trend. Indeed, Fig.^ shows 
that the slow-rotator sequence of NGG6811 tends to bend 
downwards between 0.8 and 0.9 Mq, but this trend is not 
subsequently seen in the NGG6819 periods (at the age of 
« 2.5Gyr). As this may outline some remaining problems 
with the NGG6811 data, which might result in introducing 
biases in our results, we choose not to use the NGG811 data 
in the fitting, but we verify a posteriori the compatibility 
of the final results with this data. 

The results of the fitting of the MGMG parameters are 
reported in Table[^ Fig.i and E Fig.E shows that the 
Kjjj vs. M relation obtained for the KB models produces 
the lowest linear Pearson correlation coefficient, compati¬ 
ble, within the uncertainties, with no correlation. It is, in¬ 
deed, remarkable that at f « 0.55 Gy r the slow-rot ator se¬ 
quence follows the shape predicted by Barnes (20101 for the 


^ Note that the strong ©cp mass dependence discussed in the 
following implies that, in order to constrain the TZM parameters 
at M < 0.85 Mq, data at ages 3> 0.6 Gyr would be required, 
which is not available yet. 
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M = 0.85 



M = 0.95 



M = 1.05 



Age (Gyr) 


M = 0.90 



M = 1.00 



M = 1.10 



Age (Gyr) 


Fig. 6. Comparison of the stellar envelope angular velocity evolution according to models KB (black solid line), MB (red dashed 
line), and GB (blue dot-dashed line) with observations. Mass scaling is corrected according to the MCMC ht for all models. The 
corresponding core angular velocity for each model is reported in grey with the same line style as for the envelope angular velocity. 


asymptotic slow-rotator sequence (Eq.^ with an appropri¬ 
ate rescaling of the free parameter kj.m a, solid-body rota¬ 
tion regime, this shape comparison would manifestly imply 
the mass scaling proportional to /*t, implemented in model 
KB. Although stars at 0.55Gyr have not yet reached this 
regime, our fit supports the validity of such mass scaling, at 
least in the mass range 0.85-1.1 M©. The confirmation of 
the general validity of such an empirical mass scaling awaits 
the acquisition of period data for M < 0.8 Mq at ages well 
beyond the current 0.6 Gyr limit. 


Model MB also produces a significantly lower vs. 
M correlation than the KA, KS, KK, and GB model. How¬ 
ever, the Kyj rise at M > 1.0 Mq indicates that the mass 
dependence of the MB wind braking law is less accurate 
than model KB at larger mass. 


Model KA, i.e. the original Kawaler (19881 model with 
no saturation, tends to underestimate the wind braking at 
lower mass with respect to higher mass (or vice versa). Our 
fit compensates for this behavior with a systematic vari¬ 
ation of Ky, with mass. A saturation threshold constant 
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Fig. 7. Two-zone model MCMC fitting of the slow-rotator sequence (red diamonds) on the period-color diagram (black dots for 
the slow-rotator sequence, grey for others). Parameters from model KB are used, adopting (K^) = 4.5 and correcting Tcp for likely 
biases according to the power-law fitting shown in Fig.[^ 
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Table 4. Periods derived from model KB with (Kn,) = 4.5 and 
Tcp corrected using the power-law fitting to the MCMC results 
(Fig.§ for a logarithmic-spaced age grid. 


Age 

(Myr) 

1.10 

1.05 

M/Mq 

1.00 0.95 

0.90 

0.85 

100 . 

2.65 

3.04 

4.04 

4.62 

5.22 

6.38 

125. 

2.78 

3.18 

4.18 

4.76 

5.36 

6.51 

160. 

2.95 

3.37 

4.41 

4.98 

5.55 

6.68 

200 . 

3.17 

3.59 

4.67 

5.23 

5.77 

6.87 

250. 

3.44 

3.90 

5.01 

5.55 

6.06 

7.13 

320. 

3.85 

4.33 

5.49 

6.02 

6.49 

7.51 

400. 

4.33 

4.85 

6.08 

6.59 

6.99 

7.95 

500. 

4.97 

5.52 

6.83 

7.32 

7.64 

8.54 

630. 

5.80 

6.41 

7.82 

8.30 

8.54 

9.33 

800. 

6.91 

7.57 

9.15 

9.60 

9.77 

10.42 

1000 . 

8.17 

8.90 

10.67 

11.15 

11.24 

11.75 

1250. 

9.67 

10.48 

12.49 

13.03 

13.12 

13.47 

1600. 

11.62 

12.52 

14.86 

15.54 

15.66 

15.90 

2000 . 

13.64 

14.63 

17.34 

18.19 

18.45 

18.64 

2500. 

15.93 

17.03 

20.14 

21.23 

21.68 

21.92 

3200. 

18.85 

20.06 

23.68 

25.05 

25.78 

26.20 

4000. 

21.98 

23.21 

27.32 

28.96 

29.98 

30.66 

5000. 

26.17 

26.96 

31.50 

33.40 

34.72 

35.70 


in mass, as in model KS, helps to correct this behavior, 
but only by a negligible amount for M > 1 Mq . A mass- 
dependent saturation threshold, as in Eq. (171, is somewhat 
more effective, but a clear ATu, vs. M correlation remains. 

The fitted val ues of Ki for mode l GB a re very close to 
the one given in Gallet fc Bouvier (|2013l, wh ich derives 
from the original simulation of Matt et al.|( 2012|| , but show 
a clear correlation with M. A strong correlation of the fitted 
Ki with mass (usi ng percentiles of the who le period sample) 
was also found by Gallet & Bouvier (2015|), who interpreted 
it as due to a change in magnetic configuration with mass. 

A comparison of the wind braking mass dependence of 
models KA, KB, and MB in the range 0.6 < M/Mq < 1.1 is 
shown in Fig.|^ Gomparing the latter with Fig.[^ we deduce 
that the slope of the mass dependence of both models KB 
and MB is compatible with the observations in the range 
0.85 < M/Mq < 1.0. For M > 1.0Mq, however, the slope 
of model MB is too steep and the observations favor a slope 
more similar to that of model KB. Note, furthermore, that 
the effective mass dependence of model MB depends on the 
saturation regime, although this has only a minor impact 
on the Kw fit, as in model KK. 

The resulting Tcp is largely independent of the model 
adopted (Fig.l^. The largest difference is between model 
GB, which is me only one with a wind braking law not 
strictly proportional to and all the others. Note that 
our Tcp derived from model GB is very close to that derived 


by Gallet & Bouvier] (2015) for the 90th percentile of the 
whole period distribution. As discussed below, a different 
dependence results in a different angular momentum stor¬ 
age in the radiative interior, and therefore a different cou¬ 
pling timescale, but Tcp remains a steep function of mass in 
all cases. Adopting model KB with an average {K^) = 4.5, 
we obtain a Tcp vs. M relationship which is well-represented 
by the power-law scaling: 


' cp 


^- 7 . 28 ^ 


(18) 
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In Fig.^ this relationships is compared with the Tcp ob¬ 
tained from each of the models tested in the MGMG fit. 

The evolution of the envelope and core angular velocities 
is plotted in Fig.|^ For clarity, only the results of models 
KB, MB, and GB are shown. Since the mass dependence 
has been corrected by the Kyj fit, all models with a wind 
braking law proportional to produce essentially the same 
rig and flc evolution. 

In the 0.1-2.5 Gyr age range, the fie evolution deduced 
from model GB also overlaps with that deduced from all 
the other models. Differences with the other models become 
significant after 2.5 Gyr, but remain small even at the solar 
age (see Fig.[^. Model GB model predicts, however, a larger 
angular momentum storage in the core at early ages with 
respect to the other models. 

The comparison of Figs.[2and|^explains the behavior of 
the slow-rotator sequence in the PjVi vs. (B — V) diagram. 
For M < 1.0 Mq at 0.1 Gyr, fig is well below the 


flg(t) = flohs{to)\l — 


(19) 


relationship with tg = 0-55 Gyr, and Ho bs th e corresponding 
value from the non parametric fit (Sect. 3.2 1 . The evolution 
proceeds by approaching such relationship unt il t = 0.55 
Gyr, then fig decreases more steeply than Eq. (19). In the 
time-frames of Fig.[2 this corresponds to a slow-rotator 
sequence approaching the M37 sequence between 0.1 and 
0.55 Gyr, then departing from it between 0.55 and 2.5 Gyr 
towards increasing P/y/t. The behaviour for M > IMq 
is qualitatively similar, but, in this case, fig remains quite 
close to Eq. (19) until 0.55 Gyr, then decreases significantly 
faster than Eq. (|l9[) afterwards. This explains why the slow- 
rotator sequenceror M > IMq remains close to the M37 
sequence in the P/^/i vs. {B — V) diagram until 0.55 Gyr, 
and then departs from it again from 0.55 to 2.5 Gyr toward 
higher P/y/l values. 

The comparison between fig and fig shows that a quasi¬ 
solid body rotation regime, and therefore a regime in which 
the wind braking dominates the rotational evolution, is 
achieved well after1 Gyr. Furthermore, the comparison be¬ 
tween rig and Eq. (19) with tg = 2.5 Gyr shows that a quasi- 


Skumanich evolution is achieved well after 1 Gyr. This im¬ 
plies that in order to verify a departure from the pro¬ 
portionality of the wind braking law, more data at ages well 
above 1 Gyr are needed. 

Period isochrones in the period-color diagrams are cal¬ 
culated with model KB adopting {Kyj) = 4.5 and correct¬ 
ing Tcp for likely biases according to the power-law fitting 
shown in Fig.[^ These isochrones are essentially equivalent 
to those computed using any of the wind braking models 
tested after correcting the mass dependence through the 
fitted Kyj and correcting Tcp in the same way. The compar¬ 
ison with observations, shown in Fig.[^ outlines the obser¬ 
vational origin of the small discrepancies found at the ages 
of the Pleiades and M35 for M > 1.0 Mq, where photom¬ 
etry and reddening uncertainties, particularly crucial for a 
cluster like M35, affect more significantly the results. Fig.[^ 
also shows that our model reproduces satisfactorily the pe¬ 
riod distribution in NGG6811 as well, even if it had been 
originally excluded from the MGMG fit. Even the recon¬ 
structed period at M = 0.85 Mq, where the distribution 
apparently bends towards lower periods, lies on the upper 
envelope of the observed distribution and it is therefore 
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Fig. 8. Period evolutionary tracks calculated with model KB 
adopting = 4.5 and correcting r^p for likely biases accord¬ 
ing to the power-law fitting shown in (Fig.[^. 


compatible, in a statistical sense, with the observed distri¬ 
bution. Overall, however, the evolution of the slow-rotator 
sequence is reproduced with unprecedented accuracy. 

Period isochrones recomputed for a logarithmically 
spaced age grid are reported in Table|^ and Fig.|^ Inter¬ 
polation of Table]^ provides a tool for inferring stellar ages 
of solar-like main-sequence stars from their mass and rota¬ 
tional period, effectively representing an alternative gyro- 
chronology relationship that takes the physics of the two- 
zone model for the stellar angular momentum evolution into 
account. 


6. Conclusions 


Using a two-zone model MCMC fitting to the rotational 
periods vs color data from selected open clusters, we have 
compare d new proposals for the wind brak ing law with the 
classic al Kawaler (19881 one modifi ed by Chaboyer et al. 
(19951. In our tests, we coupled the Kawaler ( 1988p brak¬ 
ing law with different hy potheses regarding the satur ation 
regime, including that of]Krishnamurthi et al. (T997|. The 
wind br a king l aws tested also included that of Gallet fc 
Bouvier (|2013 ), one derived from the work of [Barnes "fc 


IT 


and |B 


arnes 


Matt et al.)(201^ 


(20101, and the new proposal by 


The MCMC fitting has been performed on the slow- 
rotator sequence between 0.1 and 2.5 Cyr. We proposed a 
new quantitative criterion for identifying the slow-rotator 
sequence based on the symmetry of the data distribution 
around the density peak in the period-color diagram. Fol¬ 
lowing this criterion, it is possible to assign a width of the 
period distribution around the density peak using the stan¬ 
dard deviation. The value of the peak density as a function 
of colors (mass) is estimated using a non-parametric fit and, 
together with the standard deviation, constitute the data 
used in the MCMC fitting and in the definition of the like¬ 
lihood function. 

The main idea behind this work was to see how accu¬ 
rately a two-zone model with parameters constant in time 
can describe the rotational evolution on the slow-rotator se¬ 
quence. In fact, it is evident that the slow-rotator sequence 


evolves in a regular and predictable way and it is then nat¬ 
ural to explore the extent to which current models can re¬ 
produce this sort evolution. 

We find that the difficulties in reconciling the observed 
behavior of the slow-rotator sequence evolu tion between 0.1 
and 0.6 Cyr with the original Barnes ( 2003|) idea of a factor- 
ization of time and mass dependence (Meibom et al.|[2009 


2011b I can be simply attributed to a transfer of angular mo¬ 


mentum from the radiative core to the convective envelope 
in the early slow-rotator sequence evolution. The two-zone 
model MCMC fitting leads to a robust (largely indepen¬ 
dent of the wind braking model adopted) estimate of the 
core-envelope coupling timescale Tcp as a function of mass 
in the range 0.85 < M/Mq < 1.10. In this mass range we 
find that Tcp scales as on the slow-rotator sequence. 

The importance of the core angular momentum stor¬ 
age and transfer to the convective envelope, in which en¬ 
hanced viscosity that is due to (magneto-) hydrodynamical 
instabilities is expected to play an important role, has been 
outlined by many authors (see, e.g. references in Sect. 4.1). 
In this work we derive, for the first time, the empirical 
mass dependence of the core-envelope coupling timescale 
in the slow-rotator regime. Our empirical mass dependence 
of Tcp represents a useful term of comparison for the theo¬ 
retical description of angular momentum transport in stel¬ 
lar interiors. By focusing on the slow-rotator sequence, in 
which the TZM with constant parameters, as described in 
Sect. |4.1| is expected to capture the essential physics of the 
angular momentum evolution, we avoid the complications 
that arise in dealing with the fast-rotator regime, in which 
the TZM wi th constant parameters may break down (e.g. 
Brown|[2014 |. The validity of the TZM for the slow-rotator 
regime is confirmed by the successful fit to the cluster data 
currently available. The wind braking law, for which the 
empirical dependence on stellar mass and angular velocity 
still gives better resu lts than purely theoretical ones (see 
e.g. Matt et al.|2015 1, remains one of the critical aspects of 
this modeling. 

Our results outline and explain the invalidity of a fac¬ 
torization like the one expressed by Eq. (Q, on which most 
of the proposed gyro-chronology relationships are based, at 
least for the early (mass-dependent) evolution of the slow- 
rotator sequence. On the other hand, our results mean that 
the search for gyro- chronology rela tionships like the one re¬ 
cently proposed by Kovacs (2015) , which is not based on 
known physics, are not justified. Kovacsj ([2015 ) proposal 
was motivated in fact by the disagreement between gyro- 
chronology relationships based on Eq. ([^ and stellar evolu¬ 
tion isochrone fitting, but the ages assumed for the clusters 
included in our analysis, which d o not differ sig nificantly 
from the isochrone ages assumed in Kovacs (20151, are fully 
compatible with our modeling of the rotational evolution of 
the slow-rotator sequence. 

With a suitable choice of the free parameters, we tested 
the mass dependence of the models under scrutiny. We 
found that the mass scal i ng of the wind braking law impli ed 
by the models of Barnes (2010|) and Barnes & Kim (2010) is 
the most consistent with the data, supporting the proposal 
that the convective turnover timescale r is a key parameter 
in such a mass dependence. A definitive confirmation of the 
general validity of such mass dependence, however, would 
require period data of older cluster and lower stellar mass 
than is a vailable to date. T he mass scaling of the recent 
model by Matt et al. (2015) is found consistent with the 
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data for 0.85 < M < 1.00 Mq, but it decreases too steeply 
with mass for M > 1.00 Mq. All the other models consid¬ 
ered show a clear residual correlation with mass, indicating 
an incorrect mass scaling. 

Small deviations from the dependence of the wind 
brak i ng law , such as those predicted by the [Gallet fc Bou-| 
|vier| (20131 model, are also compatible with observations 
in the 0.1-2.5 Gyr range. In fact, quasi-solid body rotation, 
the regime in which the wind braking dominates the angular 
momentum evolution, is achieved after 1-2 Gyr, depending 
on mass, so that small deviations from the proportion¬ 
ality can be compensated for by a different angular mo¬ 
mentum storage in the stellar core to produce, essentially, 
the same fie evolution. After correcting the m ass scaling, 
differences between the Gallet & Bouvier (20131 model and 
models based on the fr proportionality becomes not neg¬ 
ligible after 2.5 Gyr, although they remain small even at 
the age of the Sun. A corollary of this is that verifying 
small deviations from the proportionality would require 
data for clusters much older than 2.5 Gyr and that, in the 
2-5 Gyr age range, deviations from the Skumanich law for 
M « 1 Mq are expected to be small. 

From our MGMG model fitting to the data, we derive 
period isochrones and tracks from 0.1 to 5.0 Gyr valid for 
stars with 0.85 < M < 1.10Mq. These are essentially in¬ 
dependent of any assumption regarding the mass scaling of 
the wind braking law and compatible with small deviations 
from the proportionality. By interpolating between the 
grid point, these tracks can be used to infer stellar ages 
from their mass and rotational period, providing alterna¬ 
tive gyro-chronology relationships that take the physics of 
the two-zone model into account. 
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